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Abstract: We describe the first convergent numerical method to determine static black 
hole solutions (with horizon) in 5d compactified spacetime. We obtain a family of 
solutions parametrized by the ratio of the black hole size and the size of the compact 
extra dimension. The solutions satisfy the demanding integrated first law. For small 
black holes our solutions approach the 5d Schwarzschild solution and agree very well with 
new theoretical predictions for the small corrections to thermodynamics and geometry. 
The existence of such black holes is thus established. We report on thermodynamical 
(temperature, entropy, mass and tension along the compact dimension) and geometrical 
measurements. Most interestingly, for large masses (close to the Gregory-Laflamme critical 
mass) the scheme destabilizes. We interpret this as evidence for an approach to a physical 
tachyonic instability. Using extrapolation we speculate that the system undergoes a first 
order phase transition. 
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1. Introduction 

In backgrounds with additional compact dimensions there may exist several phases of black 
objects including black-holes and black-strings. The phase transition between these phases 
raises puzzles and touches fundamental issues such as topology change, uniqueness and 
Cosmic Censorship. 

Consider for concreteness a background with a single compact dimension - M'^^^'^ x S^. 
We denote the coordinate along the compact dimension by z and the period by L. The 
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problem is characterized by a single dimensionless parameter^, e.g. the dimensionless mass, 
IX = Gi^M/ L'^~^ where Gat is the d dimensional Newton constant and M is the (asymptotic) 
mass. Gregory and Laflamme (GL) [Q, |2| discovered that a uniform black string - the d — \ 
Schwarzschild solution times a line - becomes classically unstable below a certain critical 
value ^GL- They interpreted this instability as a decay of the string to a single localized 
black hole. Their discovery has initiated intensive research^ ||3|, ^, ^, |6|, ^, ^ [l^, 11, 12 



[l3| , pl| , |T^, 1£, 17| that attempted to trace out the fate of the unstable GL string: whether 
it settles at another intermediate stable phase as advocated in ^, or whether it really 
decays to a single black hole. By now there is mounting direct evidence against the former 
possibility |^, together with additional circumstantial evidence and |15] which we 

also regard as evidence against stable non-uniform black string phase^. 

Here motivated by Q we take another route, namely we address the question: what 
happens to a small localized black hole as its mass increases (by e.g. absorption of an 
interstellar dust)? Such a black hole grows and naively one expects that there is a moment 
when its "north" and "south" poles touch. Whether this is the case or not is yet to be 
established, but it is clear that some sort of instability will show up when the poles are 
getting closer. Put differently: is there a maximal mass, beyond which the black hole 
"does not fit into the circle" and there are no stable black holes. This maximal mass 
would be analogous to the GL critical mass, and would correspond to a perturbative, 
tachyonic instability . Yet another kind of instability may occur before that maximal mass 
is reached. Once the entropy of a black hole equals the entropy of uniform black string with 
the same mass, a transition between both phases will be allowed by quantum tunneling, 
or by thermal fluctuations. This first order phase transition is slower than the classical 
perturbative instability due to tunneling suppression. 

No analytical solution for a black hole is known. Moreover, even though one can expect 
approximate analytic solutions to exist for very small black holes, the phase transition 
physics happens when the size of the black hole is comparable to the size of the compact 
dimension. Hence, in this work we take the numerical avenue. 

In our first article we considered the theoretical background for the static d- 
dimensional quasi-spherical black holes (BHs). There we outlined the goals of the numerical 
study. Prime among these goals is to establish the very existence of the static black hole 
solutions. To our knowledge, there is no direct evidence in the literature that such BHs 



do exisf^ though there are positive indications for that [14|. Among other goals is the 
study of such BH solutions in various regimes and dimensions. The ultimate and the most 
interesting aim is of course to determine the point of phase transition. 



^ Later we will use another parameter x defined in ([2.^). 
^Related research includes (l8|, H. 

^Note however that the authors of did not interpret their results either as supporting or as countering 
the conjecture. 

^Arguments like: "in the limit when the radius of a BH is small compared to the compactification radius 
the equivalence principle implies that the black hole must be similar to the 5d Schwarzschild solution" , while 
intuitive are not rigorously sufficient. In fact, this argument fails in 4d with one of the space-like directions 
being curled to a circle, as there is no stable configuration of a periodic array of point-like sources. 
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Based on recent progress 21 1 we develop a numerical scheme that allows us to find 
static axisymmetric BH solutions. Our scheme is dimension independent, provided that 
d > 4. As a first step, in this paper, we apply it to the 5d case, which is the example 
with the lowest dimension^ among spacetimes with extra dimensions. In 5d we construct 
numerically a family of static BHs, parametrized by x, which is the ratio of the size of 
the black hole to the size of the compact dimension, see (p.4|). For small values of x the 
horizon region of our solutions approaches the 5d Schwarzschild solution which can be 
considered as the 'zeroth order in a perturbative expansion' in powers of x. Moreover, in 
this limit our solutions satisfy the theoretical expectations for some next order corrections 
in this perturbative analysis thereby allowing a confirmation of a new theoretical 
method through numerical "experiment" . This establishes the existence of static higher- 
dimensional BHs and shows that the Schwarzschild solution is the smooth limit of these 
solutions for x — > 0. 

We succeed to control the accuracy of our solutions up to x < xi ~ 0.20 (corresponding 
to fii ~ 0.047). Above this limit, up to the last value X2 — 0.25 (corresponding to /i2 — 
0.074), for which our solutions do not diverge the convergence rate was very slow and the 
numerical errors were not small. These values of fi should be compared with the critical 
GL mass /Ugl — 0.070. The slowdown of convergence and eventual divergence is mainly 
seen on one of our metric functions. By examining the equations of motion for our system 
we observe a "wrong" sign in one of the equations (just like the plus sign in the following 
harmonic oscillator equation, ip" + Ld'^ip = 0), which is an indication of the presence of the 
tachyon. One could expect that the tachyonic behavior is suppressed for small x values and 
it is manifest for large x values, for which there are no static BH solutions^ ^. However, 
the tachyonic behavior influences the numerics even before that critical x value and slows 
down the convergence. We believe that the problematic variable is coupled to the tachyonic 
mode, and hence when the latter drives the former to behave pathologically, is an indication 
that the system is close to the phase transition point. 

In [^3[ we derived the d-dimensional Smarr formula, also known as the integrated first 
law, for the geometry under study (see also ||l^). It is a relation between thermodynamic 
quantities at the horizon and those at infinity, relying on the generalized Stokes formula 
and the validity of the equations of motion in the interior. This naturally suggests to use 
this formula to estimate the "overall numerical error" in our numerical implementations. 
This method comes in addition to the standard numerical tests such as convergence rate, 
constraint violation etc. While it is possible that this is a lucky coincidence (though we 
believe it is not), for our solutions the Smarr formula is satisfied with 3 — 4% accuracy. 
Moreover, it is intriguing enough that the Smarr formula is satisfied with the same 4% 
accuracy even for the problematic solutions for x > 0.20. This has to do with the fact that 

^This is maybe the lowest dimensional example, but because of a very slow asymptotic decay, it is 
certainly not the simplest to solve numerically [^. We discuss this in detail later on. 

^Consider a tachyon in a box: the mode can materialize only if its inverse mass is not less then the 
dimension of the box otherwise the mode is suppressed. 

^This is a classical "revolutionary situation" : a "poor" tachyon is suppressed until the black hole becomes 
too fat. Then the tachyon rises, gets strong and destroys the black hole, heading to a new future (to another 
phase). 
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this formula relates only 3 variables of the 4 thermodynamical variables, characterizing the 
system. It turns out that the 4th variable is somewhat decoupled from the other three. 
However, as the inaccuracy in determining this 4th variable grows with x, this slows down 
and ultimately ruins the convergence. We believe that this is the variable which is coupled 
to the tachyonic mode. 

Even though the Smarr formula is satisfied to a good accuracy, we do have some 
larger inaccuracies in the solutions. One of the fields suffers from a certain convergence 
problems, and its asymptotic behavior departs by some 30% from small x predictions. This 
is exactly the "fourth" asymptotic charge that does not appear in Smarr's formula. Due to 
its approximate decoupling it is plausible that indeed we have good accuracy for all other 
measurements. Even better, we have indications that this field is reliable for a sub-range 
of x: 0.08 < X < 0.15. 

One can question what information can be extracted from knowledge of only three 
parameters for the entire sequence of BHs (x < 0.25), and knowledge of the fourth one for 
a smaller range (0.08 ^ x < 0.15). In particular, we show that the last black hole that we 
find (at x ~ 0.25) deviates only slightly from being spherical, and moreover, its poles are 
quite distant from each other. 

In addition, one can ask whether there is a first order phase transition. We cannot 
establish this with certainty, since the entropy of our last black hole (at X2 — 0.25) is still 
larger than the entropy of the corresponding uniform black string. A naive extrapolation 
of our data to larger values of x indicates that the entropies will become equal just above 
the maximal BH that we find, namely at X3 ~ 0.26 that corresponds to /U3 ~ 0.082. It 
is rather suggestive that /iGL,/U2 and fi-s are all very close each to another. Since, all the 
numbers in the system are expected to be of the same order, this fact may be regarded as 
an indication that we have found a real phase transition. Note that since ^2 — /^gl we 
come close to a first demonstration of a failure of higher dimensional uniqueness with two 
stable phases^. Finally we note, that generally in a first order phase transition one expects 
/^GL ^ ^3 < /^2- This remains to be tested numerically. 

While we expect that the instability we found corresponds to a physical one we stress 
that we cannot rule out the conservative possibility that it is a manifestation of imperfec- 
tions of the numerics. Since our numerical scheme is independent of the dimensionality of 
the problem provided d > 4, the immediate aim for the future work would be its applica- 
tion to higher dimensions, d > 6, where the asymptotic fall off is faster and the solutions 
might be more stable^. 

In section ^ we describe our system. We employ the "conformal ansatz" and derive 
the equations of motion and the boundary conditions. A short excursion into theoretical 
background (summarized from |2^) is made in section ^. Our numerical implementation is 
described in detail in section |^, where we also describe various tests. The results are listed 
in section ^. We outline future directions in the final section |6[ In appendix we derive 

^Although we did not demonstrate that we assume that our BHs solutions are stable. 
^In fact the preliminary results show that the picture that we find in 5d is qualitatively unchanged for 
d > 6. 
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the d-dimensional field equations and boundary conditions for tlie cylinder M"^^^'^ x S^. In 
appendix ^ we consider the asymptotic behavior of the equations. 

We also refer the reader to independent work by Kudoh and Wiseman who performed 



recently related calculations in 6d |25]. 
2. Formulation 

In this section we focus on the five dimensional case ~ we derive the field equations and 
discuss the boundary conditions (b.c). Equations and b.c. on a general d-cylinder are 
discussed in appendix |^ . The fifth spatial direction is denoted by z and it is compact 
with a period L, i.e. z and 2: + L are identified . We consider static localized BHs with an 
S'^ horizon topology. We assume spherical symmetry (50(3) isometry) of the 3 extended 
spatial dimensions and we denote the 4d radial coordinate by r. 

2.1 Choice of coordinates 

We consider a static axisymmetric metric which is built out of three functions. We adopt 
a conformal (in the {r, z} plane) ansatz of the form 

ds^ = -A^dt^ + e^^{dr'^ + dz^) + e^^r^d^l , (2.1) 

where A, B and C are functions of r, z only and = dO^ + sin^ 9 dcfP' 

To describe a BH it is convenient to transform to polar coordinates, defined by 

r = />sinx, z = pcosx, (2.2) 

since the BH horizon is represented by a closed curve in the {r, z} plane. The metric in 
these coordinates reads now 

ds^ = -A^dt^ + e^^{dp^ + p^dx^) + e^^ W xd^l , (2.3) 

To simplify the numerical procedure it is desirable that the boundaries of the inte- 
gration domain^'' lie along the coordinate lines. Note that by choosing the ansatz (p.lj), 
or ( |2.3D we still did not fix the gauge completely. There is still a freedom to move the 
boundaries of the integration domain by a conformal transformation. It was shown in 
1 21 that using this conformal freedom the horizon boundary could be set at a constant 



radius ph, leaving the periodic boundaries along z = const lines. Thus the domain is 
{(r, z) :\z\ < ^ z^ > Pj^}, where for future use we define the half-period L = L/2 of 

the compact circle, see Fig. (|l]). In addition, by fixing the ratio of the radius of the horizon 
to the period of the circle 

X := ^, (2.4) 

all residual gauge freedom is eliminated. In our implementation, we set ph = 1, without a 
loss of generality, and generate different solutions by varying L. 



^"What we call here "domain of integration" could be called alternatively "domain of definition", "domain 
of relaxation" etc. By this term we refer to the region of space-time where we solve our equations. 
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(a) (b) 

Figure 1: A spacelike slice of the black-hole spacetime. (a) In the {r, z} plane the black hole's 
horizon is a curve with a spherical topology. (5) There is a conformal freedom to transform the 
domain to {(r, z) : \z\ < L, + > p^}- By fixing ph/L the domain is uniquely specified |^ . 



The fact that x cannot be changed freely for a given solution implies that a; is a 
characteristic parameter analogous to the (normalized) total mass or the temperature even 
though it does not have a clear physical meaning. For example, if one enforces the horizon 
of a BH to be at a fixed radius set by some given x, it would be excessive to specify also 
the temperature. Conversely, specifying the temperature one does not have the freedom to 



constrain the location of the horizon |21|. 



In polar coordinates the reflecting boundary of the compact circle, z = 0, is at x = 7r/2, 
but the periodic boundary, z = L, does not lie along a coordinate line in the {p, x} plane. 
The treatment of this irregular boundary introduces a certain complication in the numerical 
scheme as described in section ([4.1|). Nevertheless, we believe that it is preferable to work 
in polar coordinates (|2.2| ) and to have an irregular boundary at z = L, rather than work in 
rectangular coordinates {r, z} and have an irregular boundary at the horizon. Intuitively, 
this is because we expect that the region near the horizon would become the region of the 
'activity' as x increases. 

For numerical reasons it would be convenient to use another angular coordinate 

e = cos(x) . (2.5) 

The benefit of using this coordinate is twofold. First, the irregular z = L boundary has a 
particularly simple representation, p = L/S^. Second, as we explain shortly, the coordinate 
singularity at the axis, r = 0, becomes first order instead of second order. 

2.2 Equations of motion 

Our basic equations are the five-dimensional time independent vacuum Einstein equations. 
There are five equations in 5d: two are equations of motion for A, C, while variation with 
respect to the metric in the (r, z) plane yields three additional equations. In the conformal 
ansatz one of them is an equation of motion for B while the other two result from the gauge 
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fixing. The equations can be combined in a way that three of them wih take the form of 
elhptic equations, which we cah the interior equations. The other two combinations that 
contain a hyperbohc differential operator will be termed 'the constraints'. These constraint 
equations are not independent as they are related to the interior equations via the Bianchi 
identities. 

In order to obtain the interior equations we can follow the general procedure described 
in appendix |A|, or alternatively, use a suitable symbolic math application e.g. GRTensor 
||2^ to evaluate the relevant quantities. In either route one obtains the interior equations 
which are the following combinations of the components of the Einstein tensor: Qg + 
l/2gj + l/2gP - 2gl , 2gl - 2g^ - 2gP + and g^g+g^ + g^p- gl They can be written 
respectively as 

AA + ^ {-i + (1 - e)d^C) + 2dpA Q + dpC^ = 0, (2.6) 



Here we used the variable ^ instead of x and the Laplacian becomes A = 9^ + {1/ p)dp + 



The constraint equations expand to 



^ p^ 1 A 1— ? P 



+ 2 d^B dpC - 2 d^c dpC + ^^^^^^ ^"^^ - 2 dp^c] = o (2.9) 



GP_Gi= gg-4 (e - 2 (1 - e) d^B) AUd^B-d^C) _ 
- Ap^ + p2 

2(1- e^) (2 d^B - d^C) d^C (1 - e^) dlA 
p2 ^ Ap^ 



2 i^-^ d^C + (1 - e) die) dpA (i + 2 dpB) 
^ ^2 + A + 

/I \ d'^A 

+ 2 (2 dpB - dpC) {- + dpC\ - ^ - 2 a^C = . (2.10) 

Assuming that the interior equations are satisfied, the Bianchi identities ga-p = 0, imply 
Q the following relations between the constraint equations 

d(U + a^V = , 
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-d^U = , (2.11) 



where ( = log p., and we define the rescaled constraints U = P\/—g {^Gp — /2 , V 
P^y/^Sp with g = det g^, 



A nice feature follows |^]. The constraints lA and V satisfy the Cauchy-Riemann ( 2.11| ) 



relations and hence each one of them is a solution of the Laplace equation. Hence, if one 
of the constraints is satisfied at all boundaries and the other at a single point along some 
boundary these constraints must be satisfied everywhere inside the domain. This fact will 
be referred hereafter as the "constraint rule" . In our implementation, following the choice 
in we imposed V = along all boundaries and Z// = in the asymptotic region. It is 
important to check and confirm that the constraint U and V are satisfied everywhere for 



our numerical solutions, as we describe in section 4.2. 
2.3 Boundary conditions and constraints 



The interior elliptic equations (2.6-2.8|) are subject to boundary conditions. In this section 
we describe the boundary conditions that define the problem completely. The integration 
domain is defined by {(r, z) : < z < L, + > p^^, designated by the thick dashed 
line in Fig. ||. The boundary conditions are specified on the axis, at the horizon, in the 
asymptotic region and at the reflecting and periodic boundaries z = L and z = 0. 

2.3.1 The z = and z = L boundaries. 

On the reflecting, z = 0, and the periodic, z = L, boundaries we impose 

d^'tp = 0, ip = A,B,C. (2.12) 

While at the reflecting boundary this condition is simply = 0, at the periodic boundary 
its implementation is not direct, see section [4.1.1| . 

2.3.2 The r = axis. 

Regularity of the metric on the axis (absence of a conical singularity) requires 

B = C. (2.13) 



We use this equation as a Dirichlet condition for B. Equation ( |2.7| ) is not solved at the 
axis but it is only monitored there. For A and C the boundary conditions are automatic 
- on axis the (interior) equations for these functions become flrst order in derivatives 
normal to the boundary and have precisely the form of a b.c. Namely these equations 
are already incorporate b.c. and these do not need to be additionally specified. We term 
this an "automatic boundary condition". This occurs because of our particular choice 
of the angular coordinate: we use ^ instead of x- The axial symmetry of the problem 
dictates the dr = condition for the metric functions, which translates to dy. = in 
spherical coordinates and \/l — = in our coordinates. But on axis ^ = 1 and hence 
this condition need not be imposed in our coordinates. While the coordinate singularity 
at the axis is quadratic, ~ sin(x)~'^ when using X; it becomes linear ~ (1 — ^)~"^ when 
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using ^. With this advantage there is, however, a drawback: the metric functions are 
not differentiable at ^ = 1. This requires a modification of the numerical scheme there, 
replacing the second order normal derivative of the interior equations by a first order one 
due to the considerations above as described in subsection ( |4.1| ). 

2.3.3 The horizon 

The horizon in our construction is located at ph = 1- For static solutions various notions 
of the horizon coincide - the event horizon (globally marginally trapped ), the apparent 
horizon (the outermost boundary of locally trapped surfaces) and the Killing horizon are 
all the same. The latter characterizes the horizon as a surface where 

A = 0. (2.14) 

This implies that along the horizon 

d^A = d^A = 0. (2.15) 

Even though the horizon normally is not singular (in curvatures) our equations do become 
singular there as the function A vanishes. Now we describe how the physical regularity 
of the equations at the horizon gives boundary conditions for our functions^^. Expanding 
Eqs. (^^,2.8) at the horizon we obtain the condition 



dpC = -1 . (2.16) 

We still need a condition for B. We obtain this condition from the zeroth law of the black- 
hole mechanics (or thermodynamics), namely that for static solutions the surface gravity 
must be constant along the horizon (see for example ||2^). The surface gravity along the 
horizon reads 

K = e-^dpA, (2.17) 
and the derivative of k along the horizon vanishes 

d^K ~ d^B - ^ = 0, at ;9 = Pi,. (2.18) 

The upshot is that the boundary condition for B can be obtained in one of the forms: 
either 

from ( 2.17 ) and ( 2.13| ), or by integrating ( p.l8| ) outwards from the axis along the horizon. 



In our implementation we used the former form. However we have checked that a corre- 
sponding solution obtained by using the other option differs only slightly from our original 
one. 

Note that the condition ( 2.18| ) implies that eqn. (|2.9D (or V = 0) is guaranteed along 
the horizon, and vice versa. 



'^We assume hereafter that 9p^|p^ 7^ 0, i.e. the horizon is not degenerate. 



-9- 



We can get a different condition for B as well. Examining eqn. ( ^.10[) (U = 0) one 
obtains the condition 

dpB = -1 , (2.20) 

which is necessary to ensure regularity of that equation along the horizon. 

Altogether we now have too many conditions at the horizon: four boundary condition 
(2.14, 2361 ,2.19,2.20) for the three metric functions. However, as explained in it is 
unnecessary to impose both constraints (U = and V = 0) at the same boundary, and 
actually it is necessary not to impose both in order to protect the problem from being over 
constrained. 

Out of ( p.l9 , 2^20| ) we choose to impose ( p. 19 ). The condition ( p. 20 ) is satisfied for these 
solutions. The error becomes smaller with grid refinement and reaches 2% for our finest 
grid. For the sake of completeness we also obtained solutions using (|2^ instead of (l2lg| ). 
However, these solutions do not have a manifestly constant surface gravity. The variation 
of K along the horizon is small for small x, but can reach as much as 15% for larger x values. 
The overall difference between the two solutions is maximal near the horizon, being of the 
same 15% magnitude. This difference fades off asymptotically and the constraints are still 
satisfied (with the same accuracy). We conclude that, in principle, it is possible to use the 
condition (2.20) to generate solutions, though the numerics should be refined further to 
reach an acceptable accuracy. 



2.3.4 The asymptotic boundary. 

Performing a Kaluza-Klein reduction one observes that the z-dependence of all fields is 
carried by massive KK- modes and hence fades off exponentially for large r. Thus, in the 
asymptotic region we can rewrite Eqs. (2.(:-|2.8) retaining only the r-dependence. Defining 
for convenience C = log{r), we get 



A" + A' + 2A'C' = 0, 



2A' 



B" -B' - 2C' - {C'f - — (C + 1) - 1 + e 



2B~2C 



A 



0, 



A' 



C" + 3C' + 2{C'Y + — (C + 1) + 1 - e 



2B-2C 



0. 



(2.21) 



Here the derivatives are calculated with respect to 

Asymptotic flatness at r ^ oo, requires A — 1 = B = C = 0. Linearizing the above 
equations we can solve them analytically (see appendix ^) with these boundary conditions 
obtaining 

, log^ r 



A 
B 
C 



1 
b 



- + 0{ 

r 



+ 0( 



1 



c\og{r) I 

r r' 



(2.22) 
(2.23) 
(2.24) 



where we also included the order of the corrections. Note the logarithmic term in C. This 
log-behavior is specific to 5d and indicates a very slow asymptotic fall-off. At leading order 
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the coefficients in (|2.22| ) are related by ||2^ 



a-2b + c = 0. (2.25) 

In the numerical solution one can impose the simplest Dirichlet conditions: A — 1 = 
B = C = at the asymptotic boundary. However, since in our numerical implementation 
the "infinity" boundary is located at a finite r this option appears to be too crude. One 
can improve that by going to the next order in the expansion (2.22) and using (2.25) to 
get the refined conditions 

-^{Ar) = 1 , (2.26) 
d(Br) 

C = {-l + A + 2B)log{r) . (2.28) 

Here we have rewritten the conditions in a form convenient for a numerical implementation. 

Unfortunately, we discovered that these linear conditions do not lead to a convergent 
scheme. To understand this, observe that the function with the slowest decay is C. In 5d, 
to resolve the difference between the first two terms in its asymptotic expansion with just 
10% accuracy one has to go to r ~ exp{10)ph - the log strikes hard! When the maximal r is 
not extremely large (which is the case here for practical reasons) the non-linear corrections 
appear to be important for stabilizing the scheme 

Recall the "constraint rule" which will help us to derive a more subtle b.c. for C. In 
accordance with it we choose to enforce V = along all boundaries. The rationale behind 
it is that this constraint is satisfied trivially at the axis and at the refiecting boundaries, 
asymptotically it decays exponentially fast, and only at the horizon this constrain is not 
trivial and yields ( 2.1S| ). The second constraint must not be imposed at the horizon. At 



the axis it vanishes. Along the refiecting, z = 0, and the periodic, z = L, boundaries this 
constraint does not carry any new information as it is just a linear combination of the 
interior equations. Hence, we are left with the asymptotic boundary. This boundary can 
be potentially dangerous since on one hand, (Op — decays here and on the other hand 
the measure P\/—g blows up. This competitive behavior can result in an unpredictable U. 
Thus, to guarantee that the constraint is satisfied, the natural and unique place to impose 
U is the infinity. 

The upshot is that instead of the linear condition ( p. 28 ) we compute C at the asymp- 



totic boundary using the constraint equation U = 0. By doing so we stabilize our algorithm 
and satisfy the "constraint rule" . Note that at the leading order the vanishing of U is con- 
sistent with the linear condition ( p. 28 ). 



3. Thermodynamical and Geometrical Variables 

In this section we briefly summarize the results from our previous paper [^] with a par- 
ticular focus on the 5d case. Hereafter we work in units such that G4 = 1. 
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At infinity 

As stated in |^3| (see also [^]) there are two energy- momentum charges that asymp- 
totically characterize our configuration. These are the total mass, m (or the dimensionless 
mass /i := m/L), and what we called the tension. The latter is what an observer at infinity 
interprets as the tension of an imaginary string stretched along the compact circle. These 
charges can be calculated in terms of the numerical asymptotics |23| 



L 





1 


2 


-1 




a 


T 


~ 2 


1 


-2 




b 



(3.1) 



The asymptotic mass can be calculated also in a direct way, using the free energy or the 
Hawking-Horowitz prescription^'^ |27| which gives the same result. In the linear regime. 



using the relation ( 2.25 ) in the above formulae one can express the physical charges in 
terms of any two of the numerical asymptotics. 



At the horizon 

The characteristic quantities at the horizon are the surface gravity (the temperature) 
and the area (the entropy) of the horizon. We have already defined the surface gravity 
in (2.17), and the temperature is proportional to it T = k/2it. The entropy is related 
to the surface area by the famous Bekenstein-Hawking formula Sbh = (^/^Gat). In our 
coordinates the surface area reads 

^3 = AttpI I' e^+2C0rr^^^ ^ An{2Lfx' j' e^+^cy^TT^^^ _ (3 2) 

Out of these thermodynamical variables a single dimensionless quantity can be formed 

A^''^ := A3 . (3.3) 

In addition to the 3- area it is useful to define a pair of 2-areas of horizon sections. The 
equatorial 2-area section is given by 



^11 = 47Tple^^. 

The 2-area of the section of the horizon along the axis is just 

'1 



(3.4) 



(3.5) 



1 



With these 2-areas we define the eccentricity^^ or the "deformation" of the horizon as 



^11 



(3.6) 



^^The Hawking-Horowitz mass coincides with the ADM mass when both are apphcable 

^^This definition differs from the standard definition of an ellipse's eccentricity. It is analogous to a/b—1 = 



(1 



=2^-l/2 



1 + 0.5 e , where e is the conventionally defined eccentricity. 



- 12 - 



Finally we define "the inter-polar distance" which is the proper distance between the 
"north" and the "south" poles of the black hole calculated along the axis. This distance 
reads 



-^poles ~ ^ 



dze , at r = 0. 



(3.7) 



Ph 



Small black holes 

For small black holes (x <C 1) the tension should vanish, r ~ 0, according to Myers 
1 19 ] . In this case we have 



a 




4/3 


b 




2/3 



m, 



(3.8) 



In addition, in this limit we have 

x^. (3.9) 
Small black holes are expected to resemble a 5d Schwarzschild black hole for which we have 



5d Schwarzschild : 



ASch 

^3 



Wtt-'pI Af'' = 16TTpi, K = ^. (3.10) 



More generally, the dimensionless variables can be expanded in a Taylor series as a function 
of X. It can be shown [^] that this expansion for A^'^^ takes the form 



= 27r2(l-3-C(2)x2 + ...) , 

where C,{2) = 7r^/6 is the Riemann zeta function. 

The analogous expansion for the eccentricity reads 



6 = ^C(4)x4 + 



(3.11) 



(3.12) 



Thus the prediction is that e is positive, i.e. the black hole becomes prolate along the axis. 
This agrees with the intuitive expectation that the black hole should approach a string 
shape as it grows. 

Other dimensionless quantities are the 3-area 



^3:=^=2vrV + ...^ 



the surface gravity and the temperature 
R := kL = + . 

and the polar distance 



and T = 2ttx'^ + . , 



L 



poles/-^' 



(3.13) 

(3.14) 
(3.15) 



Smarr's formula 

The Smarr formula, also known as the integrated first law, is a relation between the 
thermodynamical variables of the problem both at the horizon and at infinity. It can 



-13- 



be obtained either from the (differential) first law together with scale invariance, or by 
computing the Gibbons-Hawking free energy and combining it with the expression for the 
massi"^. We find pl 

A-i = —aL. (3.16) 

This formula relates the horizon characteristics A3 and k with the asymptotic variable a, 
together with the dimensionful parameter L. This formula is an important test for our 
numerical solutions. 



A phase transition 

One of the most important questions that we aim to answer is whether there is a maxi- 
mal (dimensionless) mass of the black hole phase as anticipated in ||8| . This is analogous to 
the minimal mass fj,GL of the uniform black string below which the string is classically un- 
stable. What happens to a black hole more massive than the critical black hole is unknown 
and constitutes one of the puzzles of this system. The appearance of a critical mass should 
be signaled in the numerics by a very slow convergence and ultimately no convergence. 

Given the asymptotic mass ( |3.lD of the black hole we can calculate the area of the 
corresponding black string of the same mass 

Abs = 4:TT{2GimfL. (3.17) 

or in a dimensionless form 

Abs ■■=^ = 16V- (3.18) 

While the existence of a maximal mass designates a perturbative (tachyonic) instability, 
the solution with Abs = A3 designates the point of the first order transition between the 
black hole and the black string phases. This transition can occur quantum mechanically, 
via tunneling, or by thermal fluctuations. 



Summary 

In our problem we define four thermodynamic measurable quantities namely the hori- 
zon 3-area, the surface gravity, the asymptotic mass, and the tension, as well as two 
geometric quantities the eccentricity of the horizon, and the inter-polar distance. The 
thermodynamic ones are related by the very non-trivial Smarr formula ( |3.16| ) . The validity 
of this formula for our measurables is one of the most important tests for the numerics (in 
addition to usual numerical tests of convergence, constraint violation etc.) This is because 
the Smarr formula relates horizon variables, with asymptotic ones. Hence the degree of 
violation of this formula can serve as an indication of the global accuracy of the numeri- 
cal method. Another important task is to check whether the measurables have the small 
X asymptotics as expected/derived theoretically. We use dimensionless measurables: in 
this form the relations between them remain the same regardless of whether or L are 
varied^^ . 

^*See for the appearance of the scalar charge in the first law 

^^For a fixed numerical lattice spacing it is not the same to vary ph or L. 
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4. Numerical Implementation 



In this section we describe our numerical algorithm for solving the system of partial non- 
linear elliptic equations ( |2.6[j2.8| ). We estimate the rate of convergence of the numerics 
and check that the numerical errors are small and the constraints are satisfied. All our 
simulations were written in Fortran. The typical run time on a 2GHz Pentium4 PC took 
about 1-2 days. The output of the code was analyzed and visualized using Matlab. 

4.1 The scheme 

The numerical technique that is often implemented to solve partial elliptic equations is an 
iterative method, called 'Relaxation', see e.g. |3^, In this method the solutions are 
iteratively corrected, starting from some initial guess, until a desired accuracy is reached. 
For non-linear equations a modification is needed. Often an iterative Newton procedure is 



combined with relaxation to find solutions of non-linear equations [pO| , 31]. 



4.1.1 Numerical lattice and discretization 

Near the horizon we employ polar coordinates {p,^}- Asymptotically, however, cylindrical 
coordinates are the natural choice. In order to use both we choose to divide our integration 
domain into two parts: (i) 'The nearby region' near the horizon is covered by polar coordi- 
nates, (ii) 'The asymptotic region' is glued to the nearby domain from the outer, far side 
and is covered by cylindrical coordinates. The two patches overlap in order to exchange 
information about the functions during the relaxation. 

We discretize our equations on a lattice that covers the domain of the integration. We 
employ finite difference approximation (FDA) in which one replaces derivative operators 
by their discrete counterparts. The discrete operators are obtained by a formal Taylor 
expansion of functions at the grid points. We use an FDA which is second order in the 
grid spacing. For example, if h is the stepsize in, say, the p direction, which is sampled by 
index j, then the first and the second derivatives of a function ip at the lattice point {k,j) 
are written to second order as 

Analogous expressions can be found for all other derivatives. 

Our second order FDA incorporates a 5 point computation molecule in the interior. 
Obviously, it would be nice to have the same feature also at the boundaries. In fact, we 
retain this feature at all boundaries but the axis through the introduction of false grid 
points. For the functions that have a Dirichlet boundary condition we solve inside the 
domain using data at the boundaries. This is the method for A and B on the horizon and 
for C asymptotically. To implement a Neumann condition we introduce false grid points 
located one stepsize outside the real boundaries. Since the normal derivative is given at the 
real boundary we define the function at the false grid points using the corresponding inner 
points. For example, at the horizon: dpC = {Ckj^+i — Ckj^-i)/2Ap = —1, where jh = 1 is 
the location of the horizon p^; at the false point {k,jh — l) we have Ckj^-i = Ckj^+i + 2Ap. 
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Figure 2: 'The nearby region' of the integration domain covered by the polar coordinates p, ^. 
The thin dashed hnes mark the location of the false grid points used for numerical implementation 
of the Neumann or mixed Neumann-Dirichlet boundary conditions. 



Now we solve the equation on the real boundary just as for an interior point, using data at 
the false grid points. This is the method for C on the horizon and for A, B and C on the 
equator. The mixed Neumann-Dirichlet conditions for A and B, that are written in the form 
( 2.26 - 2!2^ ), are imposed in the same fashion. At the axis we have "automatic conditions" 
for A and C. Since the functions are not differentiable here in the ^ direction we use one- 
sided (first-order) ^-derivatives. Here our FDA, becomes 4 point and not completely second 
order as ( |4.1| ), see also the discussion in section 2.3.2| . 

The boundary z = L is rather complicated in polar coordinates while being very simple 
in cylindrical ones. This brings us to a closer examination of the coordinate patches. 

The nearby region: {p,C} - patch. The lattice that covers this domain, see Fig.(^), 
has nodes at k = 1, K^nax, -f^max + 1, in the ^ direction, and nodes at j = 0, 1, Jmik) 
in the p direction. Here Jin(^) is the coordinate of the last point that lies within the 
boundary z = L, which is represented by the curve p = L/^. The false grid points here are 
introduced by the requirement that for each inner point there will be outer points that allow 
the implementation of the regular 5 point second order FDA scheme. This implies that at 
each k the outer points occupy Jin{k) < j < Jout{k) where Jout(^) = Jm{k + 1). During 
the relaxation we sweep the lattice from k = 1 to k = -ftTmax and from j = to j = Jin (A;). 
At a given ray k we correct the outer points at A; — 1 when we reach j > Jin (A; — 1) using 
the reflection b.c. To this end, for each outer point, O, we find the corresponding inner 
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z=L 



z=0 



Figure 3: The asymptotic region is glued to the 'nearby patch'. The two patches overlap in order 
to exchange information about the functions during relaxation. 



(4.2) 



point, T, which is its mirror reflection with respect to the boundary ((r, z) — > (r, 2 L — z)). 
The functions at the inner point are obtained by a two dimensional interpolation from the 
surrounding points and then the corresponding outer point is updated. The important 
practical note is that the inner mirror points should be calculated only once prior to 
relaxation. 

We choose fixed lattice spacings both in p and ^ directions. In the p direction the grid 
is truncated at a finite /Ocutofr- For a specific stepsize, A,^, the maximal p is 

L 

Pcutoff A^ ' 

The step is chosen such that /Ocutoff ^ L (usually we took Pcutofr ~ lOL.) Note that 

there are only two grid points on this far boundary in the x direction for Pc^^off ' 

at -ftTinax — 1 and the other is at -ftTmaxj see Fig. |2|. Here the second patch of the integration 

begins. 

The asymptotic region: {r, z} - patch. This patch begins at r = r^i^i < Pcutoff and 
extends up to r = rmax- Note that there is a 'buffer zone' where both patches overlap, see 
Fig.^. The variation of the functions in this portion of the integration domain is expected 
to be small provided that Pcutoff ^ L. Thus, the lattice covering this portion does not need 
to be very dense. The grid has a simple rectangular geometry with uniform grid spacings. 
There are two false boundaries at z = 0, z = L. At the near boundary, rmin, all functions 
have Dirichlet b.c, the values being received from the 'nearby patch'. At the far boundary, 
''maxjWe implement the mixed Dirichlet-Neumann conditions ( |2. 261 - 2^271 ) for A and B and 
evaluate C from the constraint U = 0. Since for practical reasons we were obligated to take 
finite, not too large rmax (we usually chose rmax ~ O{W00)ph) the use of the non-linear 
condition for C is essential. 

The equations are relaxed on both patches one after the other. First we sweep the 
lattice of the 'nearby region' and then the lattice of the asymptotic region. Then the 
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sweeps are iterated. The patches communicate. When sweeping the near patch we use the 
information from the far patch to supply boundary conditions at /Ocutoff- For example, see 
Fig. ^, the green point at /Ocutoflf can be obtained by e.g. a bi-linear interpolation from the 
points marked by crosses. While sweeping the far patch the boundary condition along the 
stitch at Tmin come from the near patch. For example, the yellow point at rmin is obtained 
from the points marked by stars by a bi-linear interpolation. 

To relax the equations we used a scheme which incorporates Newton iteration. To this 
end, at each grid point {k,j), any function ip = A, B,C is updated according to 

r-ik.,) = r"i^.,)-.^^t^uu, (4.3) 

where £^{k,j) is the FDA equation of motion for this ip, and w is a numerical factor. 

The basic Gauss-Seidel algorithm uses lu = 1 and it leads to a very slow convergence^^. 
One way to speed up the performance is to use Successive Over Relaxation (SOR). The 
name originates from the fact that unlike the Gauss-Seidel scheme where the functions are 
corrected by exactly what they should be from the equations, in the SOR algorithm the 
correction is larger. The over-relaxation is managed by the relaxation parameter 1 < lo < 
2. Sometimes, for non-linear equations it is better to use under-relaxation. In this case 
< < 1 and the functions are under-corrected. The case uj > 1 {uj < 1) can be imagined 
as a sort of acceleration (friction). We implemented the SOR algorithm for our problem 
and found that the convergence rate is still unsatisfactory for dense grids. 

4.1.2 Multigrid technique 

The algorithm that we found to work well is what we loosely term here a Multigrid algo- 
rithm, see e.g. [^]. In the current simplified version of this method we solve the equations 
on several successive grids with doubled density. The basic idea of the Multigrid technique 
is simple - relax perturbations of different wavelengths on suitable lattices. Clearly, relax- 
ation of a long-wave perturbation on a very fine grid will require many iterations, while if 
we first relax the perturbation on a course grid and use the solution as an input for the 
fine grid there is a chance to converge faster. Another advantage of using the Multigrid 
method is that there is a natural measure of accuracy. We can compare the solutions on 
different grids and see whether and how the differences decrease, indicating convergence 
and scaling of the truncation error. We used 4 successive grids with halved stepsizes. Our 
implementation of the multigrid method was very simple. We just improved the solution 
going in one direction, namely passing from a coarse grid to a finer one. The original 
multigrid technique |31] incorporates motion in both directions allowing to relax newly 



excited modes on suitable grids. This two-directional method is expected to be much more 
effective and fast (but difficult to program.) 

The Multigrid technique was implemented only on the 'nearby patch'. The asymptotic 
patch was chosen to have fixed grid spacings. This is because the variation of the fields 
is relatively small at large r and there is no need for very dense grids. Note also that by 



^In fact, in our case this method is so slow that we could not infer that it converges at all. 



-18- 



decreasing the grid spacings in the nearby patch /Ocutoff scales according to (|4.2D . Thus for 
each more dense grid /Ocutoff is doubled. We, however, chose to keep the truncation radius 
constant, defined by the coarsest grid. In this case when the grid spacings were halved the 
number of the points on the boundary at p = Pcutoff was doubled. 

In most of our simulations typical grid spacings in the asymptotic patch were Az ~ 
0.25 — 1.0, and Ar ~ 0.075 — 0.25. The typical grid spacings in the near patch were 
Ap ~ 0.1 — 0.25 and A,^ ~ 0.08 — 0.12 for the coarsest grid. One can estimate the size of the 
lattice taking typical L ~ lOph, /Ocutoflf , 'Tmin ~ 10-L and rmax ~ lO/Ocutoff ~ lOOL ~ 1000/9/,. 
In this case, the size, x Np, of the coarsest polar grid is 10 x 1000, while the finest grid 
is 160 X 16000. The size, A'^^ x A^^^j of the asymptotic grid would have been 10 x 1000, if 
the grid spacing in the r direction were uniform. 

4.1.3 Extracting measurables 

Once we have a solution the horizon variables such as k, A^, A2 are calculated in a straight- 
forward way from ( p.lTf ), ( ^^ ) and (|3.4^ ) respectively . In order to obtain the asymptotic 
mass and the tension (^]^) we have to expand the metric functions asymptotically. We use 
fitting functions of a suitable form to obtain those coefficients. For example to find b we 
need a quadratic fit in for B in the asymptotic region. To find a a linear fit is usually 
sufficient for a good result. For C we used a fit of the form ci log(r)/r + C2/r. However, 
we were unable to find a reliable fitting for C. We associate this with the slow logarithmic 
decay. 

4.1.4 Further developments 

There are always compromises in numerics between the computation time and the accuracy 
of the calculation. This has to do with the grid density. Large grids mean small stepsizes 
and hence better accuracy provided that the FDA is stable, i.e. converges as a stepsize 
decreases. On the other hand, even if there were unlimited memory resources to store 
large arrays, such dense grids would result in extended CPU time that would be needed to 
sweep such large lattices. We tried different tricks to find a reliable compromise for this 
'CPU-time vs. accuracy' issue. 

• A non-uniform asymptotic grid. Originally we implemented the asymptotic boundary 
conditions at rmax, that is 30-50 L, or about 300-500/9/i. However, we found that those 
values of rmax are not large enough to determine the mass with sufficient accuracy, especially 
for large values of x. Since in 5d the asymptotic fall-off is slow, ( p.22| ), implementing the 
asymptotic boundary conditions at these rmax values is not accurate enough. Naturally, 
one would like to extend the grid to the larg er rmax* One way to do so is to add more 
points to the grid. However, in order to reach rmax larger by a certain factor than the 
original rmax the number of grid points must increase by the same factor. The same is true 
for the increase in the CPU time. Since we need rmax of order exp(10)p/i, this makes the 
mere increase in number of points absolutely impractical. 

We used another technique - a non-uniform grid spacings in the r-direction. The 
stepsizes were scaled in the following fashion 

Ari+i = (1 + e) Ari , i = l,2,...,Nr (4.4) 
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Here e is a small number that we usually took as 0.01 — 0.03. In this case the coordinates 
of the mesh points in the r direction form a geometric progression and it is possible to 
reach quite large values of rmax with a relatively small number of the mesh points. The 
discretized equations are now modified and the truncation error at a grid point i scales as 
0{eAri) + 0{Arf) rather than just 0{Arf), which is the case for a uniform grid with the 
spacing Ar j . Provided e is small this modification is not that different from a uniform FDA 
and does not cause problems. In the z-directions stepsizes were left uniform and hence the 
corresponding FDA derivative operators remained unchanged. 

Keeping e in the above range, in order to reach large enough r (of order exp(10)/3/i)the 
grid turns out to be large enough to slow down the convergence notably. In practice, we 
could reach only r^ax — 500 — SOOph- The log strikes again. 

• Additional relaxation near the horizon. This is needed in order to increase the 
accuracy of the calculation of k and the area of the horizon and its sections (3.2-3l^). 
Especially, the eccentricity ( |3.6| ) is sensitive to the accuracy of the area measurement. This 
relaxation operates over a finite portion of the mesh in the vicinity of the horizon. The 
boundary conditions along the horizon, axis and the equator are the same as before but 
along the outer boundary one uses just the Dirichlet boundary conditions that come from 
the main relaxation. 

Over this region the metric functions were relaxed on two additional finer grids. One 
could suspect that the relaxation over a finite region would produce a mismatch along the 
stitch, that can be imagined as a kink or a 'ripple' in the metric fields. However, we have 
checked that this is not the case, but rather the behavior of the functions was smooth. The 
maximal change in the functions relative to the previous grids occurred over a few mesh 
points near the horizon. 

In addition, in a couple of runs we relaxed the equations over the entire integration 
domain on an additional 5th grid. When the area and the surface gravity, obtained in this 
case, were compared to the ones obtained in the relaxation over only a small region near 
the horizon, we found that the difference in both results is less then 0.1 %. The gain in 
computation time was however dramatic - hours vs. days. 

• The over- (under-) relaxation parameter iv. There is no universal algorithm to find 
the optimal to that speeds up the convergence. Except for few very simple elliptic equations 
with simple boundaries there is no analytical prescription to pick such an optimal u. Often 
empirical estimates are the only way to find it. In our case the estimation of an optimal u 
is even harder, since we need an omega for each of the three functions for both patches. 

While we cannot be confident that the lo-s we have used in our computations are 
the optimal ones, we can estimate the range of w outside which the code slowed down or 
diverged. The choice lja = 1 — 1-2, ujc = 1 — 1-1, u>b = .5 — 1. in the nearby region, 
and, iOA = f-^c = 1 ~ 1-2, ujb = -02 — 0.1 in the cylindrical patch usually gave reasonable 
convergence rate. Some of these values depend on x. The most influential w is u;^ in 
the asymptotic region: a slight deviation of its value from the narrow range ruins the 
convergence. 
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4.2 Testing the numerics 

As usual in numerics one has to convince oneself that a particular numerical method 
produces trustable results. Before discussing our findings, we present in this section the 
evidence that our method performs well. We show that the numerical errors decrease 
sufficiently fast indicating a global convergence of the scheme and that the residual errors 
of the equations are small. In addition, in GR one has to ensure that the constraint 
equations are satisfied. We show that they are satisfied to a good extent. However, this 
ceases to be the case when x is 'too large'. For x above a certain value xi ~ 0.20, the 
convergence is slowed down, the constraints are violated significantly and the errors are not 
small for the results to be reliable. Additional accuracy estimates come from the Smarr 
formula, which our solutions satisfy with very good accuracy. 

Numerical tests 

We relax the equations on four grids with increasing density. For the first and coars- 
est grid we give some initial guess while when moving to a finer grid we start with the 
solution relaxed on the previous grid. The first thing that a good method must satisfy is 
independence of the initial guess. To achieve faster convergence we regularly used as the 
initial guess the uncompactified 5d Schwarzschild solution restricted to \z\ < L. However, 
we checked that other initial guesses, such as fiat spacetime glued to the horizon etc. relax 
to the same final solution. As an indicator of the accuracy during the relaxation we use 
the accumulative residual error defined as 

Res ip := -^Zi^l'^i^kj - Srcip^j] >with ^ = A,B,C , (4.5) 

where n is the grid number (n = 1 is the coarsest) and the factor 4"~^ roughly compensates 
for the increase in the total number of grid points. The iterations on a particular grid were 
stopped when the residuals were reduced by a desired factor relative to some number, that 
we usually took as the initial residual calculated before relaxation on that grid. In Fig. 
^ we depict the residuals on each grid. Their behavior suggests convergence. Note that 
the decrease of residuals is not monotonic all the way down. This implies that there are 
modes that are continuously excited and then decay. For small values of x these modes 
are harmless, but have an imprint on the residuals' decay - the oscillations. For larger x 
values these modes are not suppressed - the oscillation increases their amplitude and then 
they finally diverge, see Fig.^. 

In addition to the convergence per grid we can use the benefits of the multigrid tech- 
nique and check the convergence when moving between grids. We examine how much the 
solution is corrected when relaxed on different grids. To this end in Fig ^ we depict the 
differences between the solution on the nth grid and the solution obtained on the n — 1, 
coarser grid. Since we used four grids there are three such pairs. We observe that the 
solution is corrected less on finer grids, as expected for a convergent method. Most of the 
corrections occur in the regions near the horizon and near the axis. In fact, this is the 
sort of behavior that allowed us to perform further relaxation with confidence even on finer 
grids in the vicinity of the horizon as described in the previous section. 
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Figure 4: A log-log plot of the normalized residuals, Res ij) vs. number of iterations, for x — 1/7, 
implying convergence. 



We can also estimate the rate of convergence. Assuming that the solution converges 




Figure 5: A log-log plot of the normalized residuals vs. number of iterations for x > 0.20. After 
an initial convergence the solution diverges. 
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Figure 6: The corrections (or differences) in the elliptic equations when relaxed on different grids 
for X = 0.1. One can see that the finer the grid, the less the solution is corrected. 



to some /* in the limit when the grid spacing goes to zero, we can write for the nth grid 

r = U + 0{K), (4.6) 

where /i„ is the grid-spacing for this mesh. The rate of convergence is defined by the power 
p. If p = 1 one speaks about linear convergence, \i p = 2 there is a quadratic convergence 
and so on. Taking differences between solutions on different grids and considering the 
fraction of these differences the value of p can be estimated. Of course, the value oi p may 
vary at different points of the numerical lattice. However, one can calculate the minimal 
convergence rate by taking the minimum of those p. We find that in the asymptotic region 
the minimal convergence rates are p^i ~ 3, pB ^ 2, pc ^ 'i for the three metric functions. 
In the nearby region the convergence rate is also found to be at least quadratic for all 
functions. In the above estimates of p we used only the three finest grids. The reason not 
to include the first grid is simple: its prime role is to perform a rough adjustment of the 
initial guess to the given boundary conditions. Hence, one expects that the solution on 
this grid is only a very crude approximation to the final solution. To guide the reader we 
plot in Fig. the metric functions in the asymptotic region at the equator for four grids. 
The nice convergence there can be easily seen. When this is the picture we infer that our 
method converges nicely. 
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Figure 7: A run with x — 0.12. Values of the functions at the equator in the asymptotic region 
for 4 grids. There is clear convergence. 




Figure 8: Four measurables and two asymptotics as a function of grid number for x — 1/12. 
While the 3-area, k and a converge nicely for all the a;'s that we relaxed, b does not. The absolute 
variation of the variables is small however. Note, however that the asymptotic charges, fi and r, 
converge as well. 
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Convergence of this kind occurred for intermediate values of x, roughly for 0.08 ^ x < 
0.15. Outside this range the rate of convergence is still very good for all measurables but 
one. To envisage our point it is more convenient to use the numerical asymptotics a and 
b, instead of our 'physical' measurables, the asymptotic mass and tension^"^. The typical 
behavior for our measurables as a function of the grid spacing is depicted in Fig. ^. While 
A^, K and a are observed to reach their asymptotic values, b is special - it does not seem to 
converge to a definite value. Note that the asymptotic charges, and r seems to settle to 
definite values as well. Moreover, even though b does not behaves monotonically, for small 
X its value is limited to be within a narrow range, see Fig. ^. We will see below that the 
real problems of convergence begin when the fluctuations of b are not small anymore. We 
refer to these fluctuations and the lack of accuracy in the measurements of b as the "the 
b-problem" . 

To get an idea of the overall accuracy of the numerics in the final solution, i.e. after 
relaxation on all four grids, we plot in Fig. |^ the normalized error in our elliptic equations. 
This error is defined at each mesh point {k,j) as 

SAj = ^ ^ (4.7) 



+ 



+ \Src if) 



One observes that the relative errors are very small being less then 0.02%. As x approaches 
0.20 the errors grow to a level of few percent. 

Constraint equations 

Additional insight into the accuracy of the method can be gained by studying the 
behavior of the constraint equations. It is clear that these must be satisfied for the actual 



solution of the Einstein equations. In Fig. IC we plot the absolute value of the constraints 



on both grid patches. The figure shows that the constraints are not small in the far region 
of the polar patch. We believe that this loss of accuracy is due to the geometric pathology of 
the polar coordinates in the asymptotic region: the uniform grid cells in polar coordinates 
become very thin and prolonged when viewed in Cartesian coordinates. This causes a loss 
of accuracy, since asymptotically there are very few mesh points in the polar patch. Hence, 
an attempt to compute derivatives in the z-direction that are required from the physical 
point of view, gives inaccurate result. When we passed from grid to grid the constraints 
were satisfied better. When evaluating the constraints in the Cartesian patch we did not 
observe any pathologies. The constraints were small and decreased fast in the asymptotic 
region. 

One gets a better insight for the constraints' accuracy from examining the relative 



errors in them. These are defined similarly to Eqn. (4^) and plotted in Fig. 11. One 
learns that the relative errors in the polar patch are indeed very small, being less then 
one percent, even though the absolute value of the constrains is not small. Thus, both 
Figs. |l^ and ^ are complimentary and bring to light different aspects of the constraints 



'^Note that 6 does have a physical meaning - it is the scalar (or dilatonic) charge. 
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Figure 9: The normalized error in the elhptic equations, for x = 1/9. This error defined in 4.7. It 
is encouraging that the maximal error is less then 0.02 %. We plotted here errors in both patches. 



behavior. The relative errors in the Cartesian patch are difficult to calculate because of 
the small absolute values of the constraints there. The absolute errors are small there and 
vary slowly, see Fig.|l^. Hence an attempt to evaluate the relative accuracy according to 
( |4.7] ) fails, producing unpredictable results because of roundoff errors. 

The fast convergence, small errors and small constraint violation are all attributes 
of the small x runs. For certain x the convergence rate became noticeably slower and the 
errors became uncontrollable. This was the x when the "problematic" asymptotic b started 
to dominate. Even though we may expect positive tension r > 0, namely b < a/2 in 5d b 
reached a/2 and continued to grow for larger x. This designated the maximal x for which 
we could fully trust our results, which is about x ~ 0.20. 

Applications of Smarr's formula 

Aside from numerical tests the Smarr formula ( p. 16 ) provides an important theoretical 
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Figure 10: Absolute value of the constraint equations on both grid patches. While there is a loss 
of accuracy near the far boundary in the polar patch, both constraints are small in the Cartesian 
patch. 



constraint. We found that Smarr's formula is satisfied within 3 — 4%(!), see Fig. |T^. We 
find a small systematic error: when evaluating A3n/{8TraL), which should equal unity, we 
find that the mean value of the numerical points is about 0.97 with a mean spread of less 
then 2%. This shows that our numerics produces systematically under-estimated values 
with a relatively narrow spread. In addition, there is a slight increase in accuracy when 
the asymptotic boundary is moved farther away. In this case the center of the distribution 
moves toward unity, though very slowly. 

It is intriguing that the Smarr formula is satisfied with good accuracy for all our 
solutions, including those with a problematic b. Even for x > 0.20 when the "b-dominance" 



triggered convergence problems this highly nontrivial formula continued to hold, see Fig. 12 



Together with Fig.^ this suggests that despite the fact that we could not determine the 
scalar charge accurately, the other three measurables are determined with good accuracy. 
Our working assumption will be that these measurables are trustable, even for x > 0.20 
up to the last convergent solution, at X2 — 0.25, though in this regime the assumption 
becomes somewhat speculative. 
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Figure 11: The relative errors in tlie constraints in tlie polar patch. These errors are small, being 
less then 2%, even though the absolute values of the constraints plotted in Fig.nffl are not. 
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Figure 12: Testing the Smarr formula for two values of rmax- The formula predicts A3K/(87raL) = 
1. The mean value for distribution of point designated by diamonds is .967 with standard deviation 
of .017. The same for the stars is .97 and .02 correspondingly. 



Let us discuss one more aspect of the "b-problem" . As we noted before even for small x 
b did not really converge, see Fig. |^. In addition, assuming positive tension, r > 0, implies 
b < a/2, and the equality is obtained for small black holes that have vanishing tension. 
Since 6 = for a uniform string this suggests that for small x b/a should equal ~ 1/2 
and should decrease to zero as x increases. In Fig. 13 we plot the ratio b/a. One observes 
that for small x values the points are distributed around 1/3 rather then around 1/2. In 
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accordance with our working assumption a is calculated correctly, hence we estimate that 
the accuracy of the b calculation constitutes 30%, for small values of x. On the other hand, 
we have seen that for intermediate 0.08 ^ x < 0.15, b does converge well. Hence, since the 
ratio b/a is expected to decrease as x increases, the measured value of 6/a ~ 0.3 in this x 
range can be trusted. 



5. Results 



One of the most important results of this work is that our numerical solutions are the first 
strong evidence for the existence of static black holes in a non maximally symmetric higher 
dimensional spacetime. We constructed a family of solutions parametrized by x, defined 
in (^.4|). The horizon region of the solutions in this branch tend to the 5d Schwarzschild 
solution in the x — > limit and become unstable for x 0.25. Even before that, at 
X ~ 0.20, numerical errors were not small anymore. Our analysis is not capable to settle 
with certainty what has destabilized the algorithm. We tend to interpret this as originating 
from a physical tachyonic instability, though currently we are not able to determine the 
exact location of the transition point. 



The geometry of the black holes 

Let us examine the spacetime structure of a typical member of our family of solutions, 
with X ~ 1/7. In Fig. 14 we use contour plots to visualize the behavior of the metric 
functions and gain some insight for the geometry. The function A vanishes at the horizon 
and it approaches unity asymptotically. The functions B and C decay smoothly from a 
finite value at the horizon to zero at infinity. One observes as well that the z-dependence 
disappears fast as one gets away from a black hole. In addition, one learns that the contours 
intersect all the boundaries at an angle of 90 degrees: the periodic at z = L, the reflecting 
at z = 0, and the r = axis. This shows that the numerical scheme performs well at the 
boundaries. 
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Figure 13: The ratio b/a is expected to start at 0.5 for x ^ 1, and decrease for larger x. Since 
6 = for uniform black string, one might speculate that 5/a — > as x increases. Here we see a 
different behavior. This behavior is not reliable since b is problematic (at 30% level) even at small 

X. 
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Figure 14: Contours of the metric functions for x = 1/7. The black hole is located at the left 
bottom corner. The vertical periodic direction marked by z. Note the change of the topology of 
the contours. Near the horizon they are spherical while far away the contours become cylindrical, 
indicating translation invariance along z at infinity. 



The proper deformation or eccentricity of the horizon (3.6) is depicted in Fig. 15 



Below X ~ 0.05 the values that we obtain are randomly distributed around zero with 
magnitudes of 10~^, hence we concluded that the black hole is spherical to better then 
10~^ at this regime. The main tendency to be noted in the figure is that e > 0. This means 
that as X grows the black hole becomes prolate along the axis, tending to a string- like 
form, as one could expect intuitively. Another interesting feature is that the last black 
hole that we have obtained (with x ~ 0.25) is only slightly deformed, with e ~ 6.5 • 10^'^. 
Comparing with the theoretical quantitative prediction ( 3.12| ), where |C(4) — 2.89, we 
see that the numerical "experiment" agrees exactly. Moreover, the axis intersection value 
—3.5 • 10~^ agrees with the expected zero up to the numerical errors, see Table ||. We note 
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Figure 15: Eccentricity or deformation of the horizon. It stays very small up to a; ~ 0.25, where 
our code becomes unstable. The coefficient of proportionality agrees exactly with ( ^.12| )! 

that although the agreement looks "too good", a fairly good agreement persists also for 
fits on smaller neighborhoods of x = 0. 
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Figure 16: The normalized inter-polar proper distance starts at 1 for small x and decreases as x 
grows. Note the surprisingly small rate of decrease. The insert contains a speculative extrapolation. 
If the latter is correct, the "north" and the "south" poles of the black hole will touch at x ~ 0.4. 
The two extrapolation lines correspond to a spline (dash-dotted) and to a 8-degree polynomial 
(dashed) . 

The (normalized) inter-polar distance ( p.I5[ ) is plotted in Fig. One observes that i 
decreases so that the "north" and the "south" poles of the black hole approach each other. 
This decrease however is slow, such that just before the last x ~ 0.25, £ is still very far from 
vanishing. This suggests several possibilities: (1) Our numerical solution crashes due to 
uncontrollable errors, and hence at x ~ 0.25 there is not really a tachyonic instability. In 
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this case the possibihty that I wih shrink to zero and the phase transition wih be smooth 
is not precluded by our analysis. (2) The point x ~ 0.25 is in the vicinity of the real 
phase transition and i does not drop to zero, signaling a non-smooth phase transition. 
This possibility was advocated in |^]. It is interesting to note that such features as the 
finiteness of I and the smallness of the eccentricity, just before the transition, were obtained 
in 1 14] for momentary-static black hole solutions. The small x behavior of I is surprising. 



Neglecting the effect of the black hole on the spacetime metric one would expect a linear 
decrease. However, as the figure shows the decrease is much slower^^. This phenomena is 
not understood yet. It looks like as if the mass of the BH expands space in such a way 
that compensates almost exactly for its size. This effect can be called an "Archimedes law 
for caged black holes". 



In an insert in Fig. 16 we plotted a possible extrapolation of the measured i. If this 



extrapolation is correct then the poles will touch for x ~ 0.4. 

Thermodynamical variables 

Our prime thermodynamical variables are depicted in Fig. |l^. We see that in the 
small-x limit all variables tend to their Schwarzschild values, designated by the thick dashed 
line. The uncompactified Schwarzschild solution appears to be a smooth limit of the near 
horizon region of the caged black holes under discussion. One notes that three of the four 
thermodynamical variables have a smooth behavior all the way up to a; ~ 0.25, while the 
fourth one (the tension) is somewhat less consistent and for x > 0.20 its behavior is strange. 
This is in agreement with our observations in section We have argued, based on the 
success of the Smarr formula, that the three prime measurables As,T,a are robust, while 
the fourth variable, b, has an uncertainty of about 30% based on convergence problems and 
disagreement for small x. This rather large error is presumably confined to b due to an 
approximate decoupling of equations in the asymptotic region, see appendix Moreover, 
it is this decoupling that allows the success of Smarr's formula. Due to different relative 
weights that the asymptotics a and b have in the mass and the tension calculation, see 



{3A), ^ appears smooth in Fig. 17, while r does not. 

For small x <C 1 one can consider a Taylor expansion of the thermodynamical variables 
in powers of x. We can attempt to extract the expansion coefficients by fitting the mea- 
surables with polynomials. We applied a fitting procedure for the first 8 to 10 solutions for 
which X < 0.11. The expansion coefficients together with the numerical errors are given in 
Table |l|, as well as the expansion coefficients of A^'^\ defined in Eq. ( |3.3D and plotted in 
Fig.|l8[ For completeness, the corresponding data for e is also added to that table. Every- 
where the numerical coefficients are given with the estimated error. This is defined here 
by the variation of the coefficients of the fitting functions, while allowing the fit to vary 
within the 95% confidence range, which is illustrated in Fig. 18. 

We can compare now the expansion coefficients to the theoretical predications that 
are summarized in section ^. The leading expansion coefficients, which are listed in the 
last column of Table |l|, are confirmed with a good confidence, aside for the mass, for 



^It is hard to fit but seems quadric. 
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Figure 17: The dimensionless mass, tension, 3-area and ttie temperature along the found branch. 
The dashed line designates the corresponding variables for the 5d Schwarzschild solution. The 
shown equations are an approximation for small x, the coefficients with the fitting error are given 
in Table. |. 
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t igure 18: A('')/(27r2) for small x. The leading correction agrees well with ( |3.1l| ) ! The confidence 
bounds are depicted. No other Taylor coefficients could be extracted reliably. 
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which the numeric and the theoretical numbers differ by about 20%. This is the imprint of 
the poor b behavior, since 6 is a part of the formula ( |3.1[ ) for the mass. The higher order 
corrections to e and A^'^^ match beautifully with new theoretical results ||2^, Other 
coefficients have a somewhat greater uncertainty (sometimes tens of percents). Since we 
do not have a theoretical insight, and because of the large numerical errors we regard them 
as tentative. 

We conclude that in the small x regime most of our measurables are very robust. We 
argue as well, that the problematic b affects n, but it is even more influential in the tension 
calculation, see Eqn.(3T). The latter is depicted in Fig. 17. In this plot the demarcation 
line corresponds to x ~ 0.20. In our numerical relaxation we observed a slow down of 
convergence and a loss of accuracy when approaching x ~ 0.20. This fact finds its vivid 
representation in the tension plot. The points beyond the demarcation line drop suddenly 
and the tension vanishes at x ~ 0.22.^^ 

Based on the success of the integrated first law our working assumption is that the 
entropy, temperature and the numerical asymptotic a are robust not only for small x but 



for all the solutions, up to a; ~ 0.25, see the discussion in the end of previous section 4.2 



Moreover, one may observe that even though the mass evaluation is very sensitive to b the 
main contributor is still a: fj, = {a — b/2)/L. Hence, as long as b < a/2, which is the case 
here, see Fig.^, one may assume that the mass calculation is accurate within 20 — 25% 
limits. This will be our second, more speculative, assumption. Using it, one can question 
what additional information can be extracted from out data. Addressing this question we 
ask: is there 



A phase transition? 

We interpret the instability of the numerics as a manifestation of a real physical tachy- 
onic instability, which slows down our scheme for x > xi ~ 0.20, and which finally ruins 
it completely at X2 ~ 0.25. Examining the (dimensionless) masses corresponding to the 
above two x-values, we find that fii ~ 0.047 and fi2 ~ 0.074 are not that far from the GL 
critical mass, fiQL — 0.070. Since a priori, the various instability masses in the system are 
expected to be of the same order, this coincidence is rather suggestive. This is another 
evidence to our assumption that we are approaching the real physical instability. 

^^We expect that the tension, much like the mass, is always positive Q so we regard this behavior 
of T as fictions resulting from the loss of accuracy in the measurable b. 





Fitting Formula 


/l±5/l 


/2 ± 5/2 


Theoretical /i 






1.48 ± .07 


-1.3 ±0.6 


37r/8 ~ 1.18 


^3 


fix^ + f2X^ 


19.54 ± .08 


5.6 ± 1.0 


27r2 ~ 19.74 


T-i 


fix + /2X2 


6.25 ± .05 


2.0 ± .12 


27r ~ 6.28 


AW/(27r2) 


fl + f2X^ 


.99 ± .004 


-5.18 ± .4 


fl = 1, f2^ -4.93 


e 


fl + f2X^ 


-(3.3 ±4.8) • 10-5 


2.89 ± .06 


/i = 0, /2 ~ 2.89 



Table 1: The numerically computed expansion coefficients of the thermodynamical and geometric 
variables. Theoretical predictions for the leading terms are listed in the last column. 
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Figure 19: The logarithm of the dimensionless area of the black hole, A3, and black string, 
Abs^ with the same mass. The continuous lines show extrapolation of the data to larger x values, 
suggesting an intersection and a first order phase transition at a; ~ 0.26. 



The existence of a maximal mass designates perturbative, classical instability. On the 
other hand, once the entropies of the two solutions are equal for a given mass, a first order 
phase transition between the two phases can take place. This transition will occur either 



by quantum mechanical tunneling or by thermal fluctuations. In Fig. 19 we depict the 
logarithm of the dimensionless areas A, for the two phases: the black hole and the black 
string, where Abs, is computed for the same mass, fj, = n{x). Our data does not show a 
crossing of the areas. However, a naive extrapolation of the data points, marked by the 
solid lines in Fig.^, indicates an intersection just above our last BH namely at 2; ~ 0.26. 
For this value the mass is ^3 ~ 0.082, which is also of order of the critical GL mass. 

Note that the extrapolated fj.^ is slightly larger then our instability mass, /i2, while 
the opposite is expected for a first order phase transition. However, this inequality is not 
numerically significant: due to the high degree of uncertainties near the instability point 
we cannot estimate well the critical mass fi2- In addition one cannot expect that ^3 is 
evaluated accurately, since its value will depend strongly on several (not very accurate) 



last points in Fig. 19 



It is important to determine whether fj,2 < /Uql- If so there must be a third stable phase 
to which the black hole decays, for instance the stable non-uniform string However, we 
expect fi2 > fJ-GL- 



6. Future directions 

We conclude by pointing out some future directions 
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Higher dimensions d > 5. In this paper we presented a 5d numerical implementation 
of the ideas outlined in our previous paper ||2^. In principle, after some slight modi- 
fications the code can be applied to a higher dimensional problem. While we expect 
that the instability, being physical, would still be present there, we hope to improve 
the accuracy of the critical mass estimation. 

We hope that the accuracy will improve in this case due to faster asymptotic decay. 



In addition, as we described in section 4.1.4 faster asymptotic decay implies smaller 



?'max) hence smaller number of grid points and therefore faster operation of the code, 
which is now frustratingly long (1-2 days). 

Improving the performance of the method. One can try to improve further the algo- 
rithm, by e.g. using the full multigrid technique that incorporates motion up and 
down grid numbers. In order to improve the accuracy near the axis we could use 
as an angular coordinate the angle x instead of ^ = cos(x). The benefit is that the 
axis is approached faster in the x coordinates, and so there are more grid points near 
it. The drawbacks (as we explained in the text) are that the coordinate singularity 
would be second order and the periodic boundary,^ = L, would loose its particularly 
simple representation p = L/^. 

Another choice of coordinates. The metric ansatz that we used here contains three 
functions. It would be interesting to try and implement the coordinates suggested in 
and substantiated in [10, In these coordinates the number of functions reduces 



to two, and moreover they interpolate smoothly between spherical coordinates in the 
horizon region and cylindrical ones asymptotically, thus eliminating the need for the 
two coordinate patches. 
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A. Equations of motion and boundary conditions on a d-cylinder 

In this appendix we derive the equations of motion on a d-dimensional cylinder, M"^^^'^ x 
and the corresponding boundary conditions 

The most general ansatz which is time independent, time reversal symmetric ("static") 
and with axial SO{2>) symmetry is 

ds^ = - exp(2i)dt2 + dcr2(r, z) + exp(2C')dJl/_3 (A.l) 

where z and z -\- L are identified , (ir2^_3, is the d — 3 sphere. A, C are functions of (r, z) 
and da'^{r, z) is an arbitrary metric in the (r, z) plane. Note that we dropped for now the 
prefactor in front of exp(2C'). Classically the problem scales with L, and so we can set 
L = 2tt and it can always be restored by dimensional analysis. 
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We need the expression for the Ricci scalar of a fibration. For ds^ = ds^+exp(2F(x))(isy 
one has 

R = Rx + exp(-2F)i?y - dvidy + l){dFf - 2dyA(F) (A.2) 

where the Laplacian is given by A(F) = dei{g)^'^/'^ d^{g^'^ (iet{gY/'^ dpF) , grad squared is 
{dF)"^ = g^'^ d^F dyF and dx = dim(X), dy = dim(y) (we normaUzed the Ricci scalar 
hy R = Rijij so that R^d = d{d — 1)). 
The d-dimensional Ricci scalar is 

Rd = R2- 2{dAf - 2Ai + {d-3){d- 4) e'^^ 

- (d - 3) (d - 2) {dCf -2{d-3)AC-2{d-3) (dA) (dC) (A.3) 

While deriving this formula one has to be careful to consider the A, C fibration step-wise, 
and thereby get the cross term —2{d — 3){dA){dC) in addition to — 2A(A + {d — 3)C). 
The gravitational action is 

where Gn is the d-dimensional Newton constant, dV2 = \/—g(^r z) drdz, 0^-3 is the area of 



J dtdV2 e^+^'^-^^'^Rd, (A.4) 



the unit d — 3 sphere, and from now on we will drop the prefactor ^q^^Iq^ f dt. 
After integrating by parts we get 

1 = J dV2e^+'^''-^'>^[R2+{d-3){d-4)e~^^+{d-3){d-4){dC f+2{d-3){dA){dC)]. (A.5) 

Now we need to fix the metric ansatz in the 2d (r, z) space. One can use diffeomorphism 
invariance to put it in the conformal form 

ds"^ = - exp{2A)dt'^ + exp(2^)(dr2 + dz'^) + exp(2C')dJ^^ (A.6) 

The formula for the Ricci scalar of a conformally transformed metric = ex.p{2B)g°'^ 
reads |26| 

Rd = e-^^[Rd - 2{d - l)duB -{d-l){d- 2)diBdiB]. (A.7) 
Hence in 2d with the conformally flat metric we get 

R2 = -2e-^^diiB. (A.8) 
The action (without second derivatives after integration by parts) is 

1 = Jdr dz[{d - 3){d - 4)e^+2B+{d-5)C ^ e^+'^'^~^'>^ {2 OMB + 2{d - 3)9^^(7 

+ 2{d - 3)dimC + id-3)id- 4:)diCdiC)]. (A.9) 
The equations of motion are 

i: du{B + id - 3)C) + - '^^"^ - id^Cf - - - e^^-^^ = (A.IO) 

B : dii{A + {d-3)C) + {di{A + {d-3)C)f -{d-3){d-4)e^^-^^ = 0, (A.ll) 
C: du{A + B + {d- 4)C) + {ddf + {d- 4){diA){diC) 
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where the grad squared is defined here without a metric factor (diA)'^ = {diA){diA). 
We find it convenient to redefine 



e r and e 



A. 



Solving for A(yl), A{B), A(C) we obtain the equations of motion, which in {r,z} coordi- 
nates are 



AA + {d-3)— + {d- 3){drAdrC + d,Ad,C) = 0, 



AB 



{d-3){d-4) 



drCi-+ drC 

r 



id-3) 



d,C[ {d-A)d,C + 2 
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A \r 



- + drC 
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+- 



drA n 



A 



+ drC + (d - 4) 



1 - e 



2B-2C 



0, 



(A.13) 



(A.14) 



(A.15) 



where A := + d^. These equations can be transformed to polar coordinates, {p,x}j 
where we have 



AA + {d-3)^ + {d- 3)dpAdpC + 



AB 



P 

(d-3)(d-4) 



dpC 



(d - 3)d^A 



(d^c + ctgix)) = 0, 



+ a,cl +^(a^c + 2ctg(x)) 



A 



(d-4)(d-3) 1-e 



d,A ( d,C + + ^{d^C + ctg(x)) 



(A.16) 



(A.17) 
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1 
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+ {d-A) 



1 - e 



2B-2C 



p2 sin(x)2 



0. (A.18) 



Here A := + [1/ p)dp + (l//?^)^^ is the flat Laplacian in the polar coordinates. 

These are elliptic equations, and as such they are subject to boundary conditions. 
One can see that the boundary conditions are basically dimension independent. The only 
difference will appear at the boundary at infinity, because of the different fall-off rate 
of the functions. Asymptotic flatness singles out natural radial coordinate at infinity, 
r ~ rschwarzschild- In [H we found for (i > 5 



I- A 
B 
C 



+ o( 



1 



„d-4 ' 



(A.19) 
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The feature to be noted is that C is the slowest decaying function, with the same rate for 
ah dimensions d > 5. The case d = 5 is somewhat special since C decays as log(r)/r at 
the leading order. 



B. Asymptotic behavior. 



At infinity the equations become one-dimensional, as all z-dependence is washed out expo- 
nentially fast. Defining for convenience C := log(r), and retaining in ( |A.13 - AT^ ) only the 
r-dependent terms, we obtain: 

A" + {d- 4) A' + {d- 3)A'C' 



0, 



B" - B' 



{d-3){d-4) 



A' 



C (2 + C") -((i-3)— (l + C) 



(d-4)(d-3) 



(1 



^2B-2C\ 



0, 

A^ 
~A 



C" -C' + {d- 3)C" (2 + C') + ^{1 + C) + (d - 4) (1 - e2^-2^) 



0, 



1 = B = C = 0, one obtains 



where ( )' := d/d(. 

Linearizing these equations around fiat space, A 

A" + {d- 4) A' = 0, 

B" - B' - {d-3){d- 4)C' + {d-4){d- 3){B - C) - {d - 3)A' = 0, 

0. 



7, 



C" + 2{d - -)C' + 2{d -4){C -B) + A' 



(B.l) 

(B.2) 
(B.3) 



(B.4) 
(B.5) 

(B.6) 



One observes that the first equation decouples^^ from the other two and it can be solved 
to yield A ~ 1 — ae~('^~^)^ = 1 — a/r'^~^. In 5d the solutions to the other two equations 
are B = h/r and C = c^\og{r)/r together with C5 = —a + 2h [23], which can be checked 
by explicit substitution. 

For a stability analysis let us look for a solution of form B^C ~ Bo,Co ■ exp{ikr) to 
the other two equations. A tachyonic mode would be one such that Im(A;) < 0. After 
substitution in the homogeneous equations one gets the algebraic equations 

-e -ik + {d-3){d- 4) -ik{d -3){d-4) - {d- 3){d - 4) 





'Bo' 







-2{d-4) -e + 2{d - l)ik + 2{d - 4) 

which have a unique solution when the determinant of the matrix vanishes 
i (-5 + d) {-4 + d) k- {ll-7d + d^) k^ -2i {-4 + d) k^ + k^ 



0. 



0. 



The solutions of this equation are 

ki = 0, k2 = i{d — 5), 



k^ = i{d — 4), ^4 = i. 



(B.7) 

(B.8) 
(B.9) 



Since for all the modes, Im(A;) > there is no tachyon. Note that there is a massless mode 
ki = (in 5d k2 is massless as well) that corresponds to the choice C = B asymptotically. 
This massless mode can, in principle, become the unstable one due to non-linear corrections 
or numerical errors, but we have no explicit indication for this. 

^"This effective decoupling ensures that the Smarr formula is satisfied very accurately even though 6 is 
somewhat less accurate. 
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